# This script will create partial plots for the effect of each predictor variable on a BRT model
# Also, setting up to loop and batch out to the HPC
# Establish directories

in.dir = '/home1/99/jc152199/brt/data/'
setwd(in.dir)
out.dir = '/home1/99/jc152199/brt/data/'

# Load the gbm library to perform Boosted Regression Tree Analysis
necessary=c('gbm')
#check if library is installed
installed = necessary %in% installed.packages()
#if library is not installed, install it
if (length(necessary[!installed]) >=1) install.packages(necessary[!installed], dep = T)
#load the libraries
for (lib in necessary) library(lib,character.only=T)

# Read in source file for gbm.plot function NB: source file must be located in working directory

source('brt.functions.R')

# Read in data to model

model.data = read.csv('/home1/99/jc152199/MicroclimateStatisticalDownscale/ToAnalyse/MicroMacroMinMaxASCII.csv',header=T)

# Establish a list of parameters for the BRT model, including learning rate, tree complexity, and number of trees

lr.list = c(.25,.1,.05,.025,.01,.005,.001,.005)
tc.list = c(1,2,3)
tf.list = c(.1,.25,.5,.75,1)
nt.list= c(5000,10000)

# Set up a loop to run BRT model with every possible combination of above parameters

for (nt in nt.list)

	{

	for (tc in tc.list)
	
		{

		for(lr in lr.list)
	
			{
 
			# BRT from gbm package, gbm() function
			
			brt.gbm = gbm(formula = micro_max~AWAP_max+coastdist+fpcmean+fpcvar+roaddist+solar, distribution = "gaussian", data = model.data, n.trees = nt.list[1], interaction.depth = tc.list[1], shrinkage = lr.list[1], bag.fraction = 0.5, cv.folds=10)
			
			# BRT from Elith's gbm.step() function
			
			test.gbm = gbm.step(data=model.data, gbm.x= c(4:8,11), gbm.y=3, family='gaussian')
			
			# Testing gbm.plot() function from Elith
			
			gbm.plot(test.gbm,variable.no = 0, nt = brt.gbm$n.trees, plot.layout = c(3,2))
			
			}
	
		}
		
	}

#End
